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The surprising discovery of an accelerating universe led cosmol- 
ogists to posit the existence of "dark energy" — a mysterious energy 
field that permeates the universe. Understanding dark energy has be- 
come the central problem of modern cosmology. After describing the 
scientific background in depth, we formulate the task as a nonlin- 
ear inverse problem that expresses the comoving distance function in 
terms of the dark-energy equation of state. We present two classes 
of methods for making sharp statistical inferences about the equa- 
tion of state from observations of Type la Supernovae (SNe). First, 
we derive a technique for testing hypotheses about the equation of 
state that requires no assumptions about its form and can distin- 
guish among competing theories. Second, we present a framework for 
computing parametric and nonparametric estimators of the equation 
of state, with an associated assessment of uncertainty. Using our ap- 
proach, we evaluate the strength of statistical evidence for various 
competing models of dark energy. Consistent with current studies, 
we find that with the available Type la SNe data, it is not possible 
to distinguish statistically among popular dark-energy models, and 
that, in particular, there is no support in the data for rejecting a 
cosmological constant. With much more supernova data likely to be 
available in coming years (e.g., from the DOE/NASA Joint Dark En- 
ergy Mission), we address the more interesting question of whether 
future data sets will have sufficient resolution to distinguish among 
competing theories. 

1. An accelerating universe. Current models of the universe posit the 
existence of a ubiquitous energy field of unknown composition that comprises 
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about 73% of all mass-energy and yet that can only be detected through 
subtle effects. Cosmologists have dubbed this mysterious field dark energy, 
and over the past decade, it has become an accepted part of the standard 
cosmology and a focus of observational efforts. More than that: it is fair 
to say that understanding dark energy has become the central problem in 
modern cosmology. 

In the remainder of this section we explain the scientific background to the 
problem and describe the quantities that are used in what follows. In Section 
2 we discuss several techniques for making inferences about dark energy and 
examine current results. In Section 3 we describe the data we use, which 
are measurements of a particular type of exploding star. In Section 4 we 
formulate the statistical problem as a nonlinear inverse problem 



where the l^s are scalar observables, w is an unknown function, called the 
dark energy equation of state, the forward operator T is nonlinear and de- 
pends on two unknown parameters, and the e^'s are heteroskedastic, un- 
correlated noise with known variances. In Section 5 we use features of the 
forward operator T to construct hypothesis tests that can distinguish among 
competing cosmological models with minimal assumptions about w. In Sec- 
tion 6 we present a framework for computing parametric or nonparametric 
estimators of w based on equation (1) and compute resampling-based error 
bars on the estimates. In Section 7 we apply these methods to current data 
in an attempt to distinguish among competing cosmological models. Finally, 
in Section 8 we look to the future, where planned observations, both space- 
and ground-based, will produce much larger data sets that can benefit more 
fully from these techniques. 

The story of how an unobserved and unexplained source of energy came 
to be quickly and widely accepted as a dominant component of the universe 
is an essential prelude to discussing the inference problem that is the subject 
of this paper. To tell this story properly, we need to start earlier, with one 
of the observational foundations on which all current cosmological models 
rest: the universe is expanding. 

1.1. Nubble's law and the distance-redshift relation. In 1929 Edwin Hub- 
ble observed a sample of nearby galaxies and studied two quantities for each. 
The first was the galaxy's redshift.^ Light from an object that is emitted 
at one wavelength and observed at a higher (lower) wavelength is said to 
be redshifted (blueshifted). Astronomers quantify this using a dimensionless 
parameter z, called the redshift and given by 





i = 1, . . . ,n 



(2) 




z = 



'emit 



^Vesto Slipher had first measured these redshifts a decade earlier [Shpher (1917)]. 
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where Aemit is the wavelength of the hght measured in the reference frame 
of the object when the hght is emitted and Aobs is the wavelength measured 
in the reference frame of an observer at some later time. When z > 0, the 
observed light is redshifted relative to the emitted light; when z < 0, the 
observed light is blueshifted relative to the emitted light. The light from 
an object such as a galaxy contains a mixture of different wavelengths at 
different intensities, called the spectrum of the object, and redshifting (or 
blueshifting) causes the entire spectrum to be displaced by a common factor. 
Because certain features in a galaxy's spectrum, such as emission or absorp- 
tion lines, have known wavelengths and identifiable patterns, the redshift 
z can be determined very accurately from observations. Several physical 
phenomena can cause the redshift to be nonzero, but Hubble interpreted 
the redshifts as Doppler shifts caused by the relative motion of the galaxies 
with respect to Earth, ^ where z > corresponds to an object moving away 
from us and z < corresponds to an object moving toward us. (Think of an 
ambulance's siren as it passes you at high speed; the pitch is higher as it 
approaches and lower as it recedes.) 

The second quantity Hubble studied for each galaxy was its distance from 
us. This he measured himself. Determining accurate distances to faint and far 
away objects is a challenging, fundamental problem of observational astron- 
omy. The ladder of measurement methods that astronomers have devised can 
produce reasonably accurate distance estimates in overlapping ranges from 
nearby stars out to very distant galaxies. Many of these methods rely on 
having a standard candle, a class of celestial objects whose intrinsic bright- 
ness is known. If we know how bright an object really is (called its absolute 
magnitude) and how bright it appears to us (called its apparent magnitude), 
we can determine its distance from us because the intensity of a light source 
decays as the inverse square of distance. To determine the distances to the 
galaxies he observed, Hubble used as a standard candle a class of stars, 
called Cepheid variables, that fluctuate in brightness with a period that is 
a function of their average luminosities. 

Hubble combined the two measurements — redshift and distance — to pro- 
duce a plot much like that in Figure 1 [Hubble (1929)]. This shows a strong 
linear relationship 



where c is the speed of light and Hq is a fundamental parameter called the 
Hubble constant. (Note that the linear form of this relationship is in fact 



^Slipher made a similar interpretation: that the "spiral nebulae" he observed were 
moving away from the Earth. But the claim that such nebulae lay outside the Milky Way 
was still a matter of intense debate at the time. It was Hubble who settled that argument. 
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Fig. 1. A plot analogous to that published by Edwin Hubble but produced with the current 
data described in Section 3. The vertical axis shows redshift and the horizontal axis gives a 
■measure of distance. The regression line through the origin highlights the linear relationship 
in (3). 

an approximation that only holds for small redshifts and distances.) Hubble 
thus found that galaxies in all directions have redshifts that tend to increase 
in proportion to their distances, and by interpreting redshifts as Doppler 
shifts, that galaxies are moving away from us at a velocity v proportional 
to distance: 

(4) V = Hod. 

This relationship is now known as Bubble's law. A natural interpretation of 
(4) would be that Earth lies at the center of an extra-galactic explosion, but 
this would violate the principle that humans are not "privileged observers" 
of the universe [Peacock (1999)]. Hubble was left to conclude that every 
galaxy was moving away from every other and, thus, that the universe itself 
is expanding. 

If the universe is expanding, then extrapolating backward in time, the 
universe was smaller in the past, and thus more dense (with the same total 
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mass-energy packed into a smaller volume), and thus hotter. Taken back far 
enough, this would predict a point of infinite density and temperature — a 
singularity known as the "Big Bang."^ However, basic physical assumptions 
break down within about 10~^^ seconds of that singularity, so there is cur- 
rently little known about the singularity itself. What we do know is that the 
universe began small, dense, and very hot and that it cooled as it expanded. 

Extrapolating forward in time, cosmologists could then predict three pos- 
sible fates for the universe, depending on whether the gravitational pull of 
all the universe's mass-energy is sufficient to overcome the energy of the 
Big Bang that is driving the universe's expansion. In an open universe, the 
total mass-energy is insufficient, and the universe will expand forever. In a 
closed universe, the total mass-energy is sufficient to halt the expansion, and 
the universe will collapse in on itself in a "Big Crunch." On the boundary 
between these two cases is a flat universe, where gravity and the expansion 
exactly balance.^ The search to distinguish among the three cases led to more 
questions than answers (mostly due to the lack of good standard candles). A 
leading model for the universe's early expansion [Guth (1981)] engendered 
a theoretical bias favoring a flat universe; however, the mass and energy 
needed to slow the expansion could not be found. In fact, studies of galax- 
ies and galaxy clusters suggested that matter accounts for only a fraction 
of the energy density required for flatness. And observations of the Cos- 
mic Microwave Background (CMB) — relic radiation from several hundred 
thousand years after the Big Bang when the universe became transparent 
to light — ruled out a flat universe comprised entirely of matter [Smoot et 
al. (1992)]. Yet several possibilities remained. Many astronomers thought an 
open, matter-dominated universe the most likely explanation for the obser- 
vations. Others found the theoretical arguments compelling and favored a 
flat universe with either a new form of energy or a new theory of gravity to 
explain the missing energy density. 

Theoretical cosmologists had anticipated an expanding universe several 
years before Hubble analyzed his data. Alexander Friedmann and George 
Lemaitre independently used Einstein's General Relativity to derive theo- 
retical models that describe the dynamics and geometry of an expanding 
universe. These evolved into the standard model of relativistic cosmology, 



^Strictly speaking, the term Big Bang to refers to the general notion that the universe 
expanded from an early state that was enormously hotter and denser than the present 
universe. But in common parlance, it is also used to refer to the initial event, as in "years 
after the Big Bang." 

^ An analogy to an open universe is a rocket that has enough energy to escape Earth's 
gravity and fly out into space. An analogy to a closed universe is a rocket that cannot 
escape Earth's gravity and falls back to the ground (or into orbit). An analogy to a flat 
universe is a rocket with just enough energy to escape to infinity but which would not 
escape if it had any less. 
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called the Priedmann-Robert-Walker model [Peacock (1999)]. Friedmann 
introduced a dimensionless function called the scale factor, a{t), that gives 
the ratio of the size of any region at time t to the size of that region at 
the current time to [from which it follows that a(to) = !]• (Cosmologists of- 
ten use subscripts for quantities at the current age of the universe.) The 
function a{t) describes the universe's entire expansion history. Friedmann 
derived the following equation that describes how the scale factor evolves: 

(5) r44V=^P(o 



At) J 3 ' a2(t)' 

where d{t) is the time derivative of a{t), G is Newton's gravitational con- 
stant, p{t) is the total energy density at time t, and A; is a parameter describ- 
ing the curvature of spacetime, which takes values —1, 0, or 1 depending on 
whether the universe is open, flat, or closed. 

The Friedmann equation (5) connects the universe's geometry (e.g., size, 
curvature) to its content (e.g., energy density). An observer at time t who 
replicated Hubble's study would find a linear relationship between recession 
velocity and distance to nearby galaxies with slope 

a{t) 



(6) H{t) 



a{t)' 



This Hubble parameter H{t) describes the (local) relative expansion rate 
and is constant for all observers at a particular time; at the current time, 
we get the Hubble constant Hq = H{tQ). In a flat universe — when fc = in 
equation (5) — the total energy density is determined by H and is called the 
critical density pcnt'- 

If p > /Ocriti the universe is closed; if p < Pcrit) the universe is open; if p = pcriti 
the universe is flat. Cosmologists refer to these cases in terms of the scaled 
density parameter 

(8) n^-^, 

Pcrit 

giving Q> 1, < 1, and Q = 1 for closed, open, and flat respectively. At 
times, it is also useful to decompose the total energy density as a sum of con- 
tributions from various sources (e.g., p = ^matter + ^radiation H ), and in such 

cases Q has an analogous decomposition. For instance, Qm = Pmattcr/Pcrit 
gives the fractional contribution of matter to the critical density. 

Relating these theoretical expressions to Hubble's observation, the red- 
shift z can also be expressed in terms of a. As the universe expands, light 
waves traveling through it are stretched out by the expansion increasing 
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the wavelength by the factor by which the universe has expanded between 
the emission and observation times. That is, hght emitted at time temit 
with wavelength Aomit and observed at time iobs has wavelength Aobs = 
Aemit(l + z), where the redshift z is given by 

(9) ^ = §^,-'- 

"I, '■emit ) 

For example, an object observed at z = 1 emitted its light when the universe 
was half its present size. The most distant objects yet observed have z w 6, 
when the universe was one seventh its current size. This reveals the redshift 
to be cosmological in origin, the result of the universe's expansion. Because 
expansion implies a{t) > 0, redshift can be viewed as an index of time, so it 
is common for cosmologists to parameterize time-dependent functions inter- 
changeably by time or redshift, for example H[t) or H{z). 

Redshift is a measure of distance as well. Bubble's relation (4) defines a 
frame of reference in which an observer at rest in that frame sees galaxy re- 
cession velocities proportional to distance in all directions. Such an observer 
is said to be comoving with respect to the Hubble expansion. In contrast, 
an observer moving relative to this reference frame would see systematically 
higher recession velocities behind than in front. The distance measured with 
a tape-measure between two comoving observers at time t has the form 
d{t) = a{t)cr, where r is called the comoving distance between the two ob- 
servers, which we express in units of time.^ Each thus sees the other receding 
at velocity v{t) where, because r does not change with time, 

(10) v(t) = d(t) = d(t)cr = ^d{t) = H{t) d{t). 

a{t) 

This gives the distance-velocity relation, or Bubble's law. But recession ve- 
locities are not observable, redshifts are. The comoving distance between an 
observer and an object at redshift z is derived by computing tape-measure 
distances between nearby events along the line of sight to the object and ad- 
justing each such distance for the corresponding expansion of the universe. 
This gives 

dz' 



Through equation (11), Friedmann's model ties the distance-redshift rela- 
tion to the universe's geometry. By carefully measuring objects' redshifts 



Various equally valid conventions for the units of comoving distance are used in the 
literature, and the choice does not change the results. We express r in units of Hq^ , or 
time. Note that in certain units commonly used by cosmologists, time and distance have 
the same dimensions because the speed of light provides an absolute conversion between 
them. 
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and distances, it is possible to estimate the distance-redshift relation and, 
in turn, the universe's expansion history, its eventual fate, and a variety of 
fundamental cosmological parameters. This was the basic task of observa- 
tional cosmology for many years. 

But in 1998, cosmologists discovered something surprising. 

1.2. Acceleration and dark energy. In 1998 two groups of astronomers 
[Perlmutter et al. (1998); Riess et al. (1998)] estimated the distance-redshift 
relation (11) using Type la supernovae (SNe), a class of exploding stars 
whose distance can be measured with ~15% accuracy, much better than for 
other distant sources. What they found was that a{t) is increasing; in other 
words, the universe is not merely expanding, the expansion is accelerating. 

The immediate challenge for astrophysicists was verifying that the appar- 
ent acceleration is not an artifact of incorrect assumptions or misinterpre- 
tation of the data. Since the initial discovery, many more supernovae have 
been measured and with greater precision [see, e.g., Davis et al. (2007)], and 
concerns about systematic errors have been allayed, though not eliminated. 
Also, indirect supporting evidence comes from measurements of the CMB 
combined with observations of large scale structure in the distribution of 
galaxies [Boughn and Crittender (2004); Fosalba et al. (2003); Nolta et al. 
(2004); Scranton et al. (2003)] and other types of data [see, e.g., Frieman 
et al. (2008) and the references therein]. These observations all support the 
hypothesis that the universe is both flat and not comprised entirely of mat- 
ter; in fact, matter constitutes only about one quarter of the critical density. 
These results also put stronger constraints on cosmological parameters that 
in turn sharpen the results of supernova studies. Taken together, current 
data strongly rule out a nonaccelerating model in comparison to a simple 
accelerating model. 

A more fundamental challenge thus becomes explaining the acceleration.^ 
If General Relativity accurately describes physics at large scales and if the 
universe is homogenous and isotropic at large scales as commonly assumed, 
then an accelerating universe can be explained by a heretofore unknown type 
of energy acting against the pull of gravity to speed up the expansion. This 
energy is characterized by its negative pressure. In contrast, in a universe 
filled with hot gas, which has positive pressure, the energy of the gas adds to 
its gravitational account, slowing the expansion. With a negative pressure 
"fluid," the opposite occurs, causing the universe's expansion to accelerate. 
Because its source and nature are unknown and because we cannot see it 
directly, this energy field with negative pressure has been called dark energy. 



^For more on explanations of an accelerating universe, the reader is referred to Carroll 
(2003, 2001), on which much of the following discussion in this subsection is based. 
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An alternative to the existence of dark energy is that General Relativity or 
the standard cosmological models built on it do not adequately describe the 
universe at large scales. General relativity has been strongly tested within 
the solar system and nearby universe, but not on scales roughly the size of the 
current universe. Active efforts have been made to modify General Relativity 
to produce an apparent acceleration. The simplest such modifications have 
been ruled out by other data and theoretical consistency requirements. More 
viable modifications have been developed but remain controversial. 

Both lines of inquiry remain open, but for now, the evidence appears to 
strongly support both an accelerating universe and the existence of dark 
energy. This raises the question of what dark energy is. One possible an- 
swer had been introduced by Einstein long-ago for a different reason. At the 
time that Einstein had proposed his theory of General Relativity — fourteen 
years before Hubble's observations — it was widely believed that the universe 
is static. But to achieve a solution of his equations that produced a static 
(though unstable) universe, Einstein needed to introduce a "cosmological 
constant" to his equations. He had never been enthusiastic about the cos- 
mological constant because it sullied the pure beauty of his equations, and 
with Hubble's observations and subsequent evidence for an expanding uni- 
verse, Einstein ruefully withdrew the cosmological constant, allegedly calling 
it his "greatest blunder" [Gamow (1970)], possibly because he missed the 
chance to predict an expanding universe. Years later, particle physicists res- 
urrected the idea of a cosmological constant to represent the energy-density 
contribution from empty space, a so-called vacuum energy. At cosmological 
scales, their formulation is mathematically equivalent to Einstein's. 

A nonzero cosmological constant can explain an accelerating universe be- 
cause it acts against gravity and because, as the universe expands, there is 
more space and thus an increased effect. The cosmological constant model 
specifies a constant vacuum energy throughout time. This is a simple model 
that is consistent with the available data. There are two problems, however, 
that suggest a more complicated picture. First, measured values of dark en- 
ergy are smaller than theoretical predictions of vacuum energy (from quan- 
tum field theory) by 120 orders of magnitude [Frieman, Turner and Huterer 
(2008); Weinberg (1989, 2000)]. Second, according to the cosmological con- 
stant model, we live in a time when the total energy density of dark energy 
is of comparable order to the energy density of matter, which, as we will 
see below, is a surprising coincidence. Moreover, although the cosmologi- 
cal constant model has an appealing simplicity, there is no known physical 
reason to require the properties of vacuum energy to be constant in time. 
Two critical questions then are whether the available data can rule out a 
cosmological constant in favor of dynamic, time-varying, dark energy and, 
if so, what the data can tell us about how it varies. 
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Once we move beyond a cosmological constant, there are relatively few a 
priori constraints on the dark energy. We can get some information, however, 
from the fact that the universe is accelerating. Let p{t) be the universe's total 
energy density at time t. As stated above, we can decompose p{t) into a sum 
of contributions from different sources, including a dark energy contribution 
p^g. The Friedmann equation (5) can be written as 

(12) a\t) = ^a\t)p{t)-k, 

and acceleration implies that a^{t) is an increasing function. It follows that 
a?{t)p{t) must increase as well and that neither matter nor radiation can be 
responsible for the acceleration. To see the latter, notice that pmattcr oc 
because the total mass of matter in any comoving volume element remains 
constant while the element's volume increases as a . Similarly, Pradiation o^ 
because expansion redshifts light to higher wavelength and thus lower 
energy by a factor of 1/a and increases volume as a? . 

The defining feature of a cosmological constant, however, is that its en- 
ergy density remains constant, so under this model, p-^^ oc a° and Pjy^a? 
increases. More generally, we require that p^^ oc a" for n > — 2. This brings 
the coincidence problem referred to above into relief. The ratio of the matter 
density to the dark energy density is given by 

(13) Anatter ^^_(3+n)^ 

Pde 

SO matter dominated in the early universe and dark energy will dominate 
eventually. Currently, the ratio is about 1/3, which puts us at a relatively 
unusual time in the life of the universe where the two components are nearly 
balanced. Cosmologists do not like coincidences, and one hope for dynamical 
models of dark energy is that they will explain the current balance. 

1.3. The equation of state. One way to quantify the dark energy is 
through its energy density p^g. This is an intuitive quantity that plays a 
direct role in the equations of cosmological models. 

Another convenient quantity that describes dark energy is its equation of 
state, typically denoted hy w. In physics, an equation of state is a formula 
that relates several macroscopic observables of a system; an example is the 
ideal gas law, relating pressure, density, and temperature. In cosmology, dark 
energy must have a "perfect fluid" equation of state, characterized entirely 
by its energy density p and isotropic pressure p. The simplest candidate is 

(14) p^E=^PDEC^ 

where w is a, property of the dark energy field that may vary across cosmic 
history. With some abuse of terminology, the ratio w is itself called the 
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"equation of state." In general, tt; is a function, usually parameterized in 
terms of redshift as w{z), but for the cosmological constant model, w does 
not depend on z. 

For the cosmological constant model, w^z) = — 1. This model describes 
a smoothly distributed field with constant energy density and negative pres- 
sure. A slightly more general model sets w{z) = wq, where wq is a constant 
that need not equal —1. In this case, /o^j, oc a~3(i+'^o) ^nd a{t) oc t^/CsCi+i^o))^ 
so the expansion will accelerate (in a dark-energy dominated universe) if 
Wq < —1/3. U Wq = —1, the dark energy density stays constant with time; if 
u^o > — 1, it decreases; and if < — 1, it grows. Cosmologists often restrict 
the possible energy-momentum tensors with various "energy conditions" im- 
plied by candidate (and somewhat speculative) fundamental constraints on 
solutions to Einstein's equations. A commonly used such condition requires 
w > —1, and most cosmological models follow suit [Carroll et al. (2003)]. 

1.4. Measuring dark energy: Type la Supernovae. One of the challenges 
of dark energy is that it has not yet been observed directly and can only 
be measured through its subtle influence on other phenomena. Astronomers 
have developed several methods that are sensitive to the expansion history 
of the universe and thus depend, directly or indirectly, on the dark energy 
equation of state. Of these. Type la SNe currently provide the best available 
constraint on the equation of state. 

Type la SNe are thought to occur when a white dwarf star strips mass 
off an orbiting companion star until it becomes massive enough to explode. 
A white dwarf is the remnant of a low- to medium-mass star at the end 
of its life, with nuclear fusion exhausted, but as it acquires mass from its 
companion it reaches a threshold above which a supernova occurs. At its 
peak, the supernova is typically brighter than its entire host galaxy, and 
it decays in brightness over a span of days or weeks. For reasons that are 
not yet fully understood. Type la SNe all have a similar peak luminosity, 
raising the hope that they might serve as standard candles. In fact. Type 
la SNe are close to but not quite standard candles because there remains 
substantial scatter in peak luminosity among nearby SNe. But these SNe 
also exhibit a strong empirical correlation between peak luminosity and the 
time it takes them to decrease in luminosity [Phillips (1993)]. Less luminous 
SNe decay more rapidly, while more luminous SNe decay more slowly. A 
one-parameter fit reduces scatter in peak luminosity significantly. Taken to- 
gether, these features make Type la SNe valuable cosmological probes: they 
are bright enough to be detected at great distances and act as "standard- 
izable" candles for distance determination. Potential systematic errors are 
thought to be smaller than current statistical uncertainties, with the main 
sources being (i) possible intrinsic differences between Type la SNe at low 
and high redshift and (ii) uncertainty in the extinction/reddening of light 
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caused by dust [Wood-Vasey et al. (2007)]. Observations of Type la SNe at 
different redshifts can thus be used to estimate the distance-redshift rela- 
tion, and in turn the dark energy equation of state. These are the data we 
consider in this paper. 

2. Inference for dark energy. Under mild assumptions,^ we can express 
the dark energy pressure p^^ and density p^^ in terms of the co-moving 
distance r as follows: 

1 \2/ 2r"{z)' 



(15) p^^(^)=_p„j,(___) (l + (l + z) 



and 

(16) PDE(^)=Pcrit 



/ 1 



f^m(l + Z 



\Hor'{2 

where is the fractional contribution of matter to Pcriti and ' denotes 
differentiation of r with respect to z [Huterer and Turner (1999)]. Taking 
the ratio of these functions yields the so-called "reconstruction" equation 
for the equation of state w{z): 

(,^. < , _ Hl^m{l + zf + (2/3)(l + z)r"{z)/{r'{z)f ^ 

^ ' Hinm{l + zr-l/ir'{z)r 

A variety of important cosmological models can be expressed in terms 
of w, including the cosmological constant {w{z) = —1): topological defect 
models [frustrated cosmic strings u;(z) = —1/3 or domain walls w{z) = —2/3; 
Bucher and Spergel (1999)]; various quintessence models [freezing w'{z) > 0, 
Caldwell and Linder (2005), and thawing w^z) < 0, Zlatev, Wang and Stein- 
hardt (1999) and Steinhardt, Wang and Zlatev (1999)]; and even models 
which allow w{0) < —1, such as Cardassian models [Freese and Lewis (2002)] 
and phantom dark energy [Caldwell, Kamionkowski and Weinberg (2003)]. 

Several critical questions can be directly addressed with Type la SNe 
data: 

1. Are the data consistent with the cosmological constant model? 

2. If not, do the data require that the dark energy equation of state varies 
with time, and if so, how well can we estimate w? 

3. Do the data rule out any competing theoretical models? 



In this paper we follow current standard practice and assume a homogeneous, isotropic, 
and spatially flat universe where matter is nonrelativistic and where gravity is described 
by General Relativity with the Friedmann-Robertson-Walker metric. Specifically, we do 
not examine an alternative to the dark energy hypothesis, modified gravity, for which 
some of these assumptions do not hold; see, for example, Huterer and Linder (2007). For 
the remainder of this paper, we also use units where the speed of light c — 1. 
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A finding that w ^ —1 or, more generally, that w is not constant would rule 
out the simplest explanation for dark energy, vacuum energy, and would 
point the way to fundamental new physics. Eliminating some competing 
theoretical models or producing a sharp estimate of w would strongly con- 
strain the theoretical explanations for dark energy. 

The results of studies to date, combining several types of observations, 
have developed a robust consensus around several findings [Frieman, Turner 
and Huterer (2008)]. First, there is strong evidence that the universe is ac- 
celerating. Second, under the current theoretical framework, there is strong 
evidence that dark energy exists, with a critical density Oj^j, ~ 0.76. And 
finally, the best available estimates suggest that w « —1; thus, current data 
are well fit by the cosmological constant model. 

Attention is therefore focused on future Type la SNe data sets, which 
promise to be orders of magnitude larger and, consequently, to provide 
tighter constraints on w. For instance, two particular observatories will be- 
gin to collect rich samples of Type la SNe during the next decade. The Large 
Synoptic Survey Telescope (LSST)^ is a ground-based instrument that will 
scan the entire sky every few nights and is expected to detect hundreds 
of thousands of Type la SNe per year, while the DOE/NASA Joint Dark 
Energy Mission (JDEM)^'^ is a space-based instrument that will observe 
thousands of Type la SNe, many at relatively high redshift. In addition, 
additional observations of other indirect probes of dark energy may provide 
strong constraints on w that complement those produced by supernova data. 
With so much data, it is expected that the dark energy equation of state 
will be well determined. (Indeed, systematic errors may come to dominate 
statistical uncertainties in the inference problem.) But whether this is suffi- 
cient depends on the structure of w and on the subtlety of the features that 
must be determined to distinguish a theoretical explanation for dark energy. 
Statistical techniques that can efficiently capture complex structure with a 
minimum of extraneous assumptions will be needed. 

Many approaches to estimating w have been used in the literature, but 
there are three main threads. In the first thread, one assumes that w lies in 
a specific parametric family w{z;9), maps this family through the forward 
operator in equation (1) to a parametric family for the comoving distance 
r(z;9), and uses maximum likelihood or a comparable criterion to estimate 
6. Many authors utilize this approach, most recently including, for example 
[Barboza and Alcaniz (2008) and Liu et al. (2008)]. Though many parameter- 
izations have been explored, the most common include w constant, w = wq, 
and w linear in the scale parameter, w{z) = wq + wiz / {1 + z) [the so-called 
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CPL parameterization; Chevallier and Polarski (2001); Linder (2003)]. The 
number of free parameters is generally limited to two, with arguments made 
that the data will not constrain additional parameters [Linder and Huterer 
(2005)]. A variant of this approach is to fit a piecewise constant model to 
w and use the fitted covariance matrix to rotate into a basis in which the 
coefficient estimates are independent [Huterer and Cooray (2005)]. In the 
second thread, one assumes that r lies in some nonparametric class; esti- 
mates r and its derivatives by nonparametric smoothing, and uses equation 
(17) to estimate w [e.g., Daly et al. (2008)]. And in the third thread, one 
derives the forward operator in terms of the energy density instead of 
w and applies one of the foregoing procedures to estimate w (e.g., Wang 
and Mukherjee (2004)]. There have also been efforts to eschew the equa- 
tion of state representation for dark energy and directly estimate kinematic 
parameters (such as the deceleration function q, see Section 4) from which 
the expansion history can be derived [e.g., Shapiro and Turner (2006)]. Fi- 
nally, there have been several papers [e.g., Huterer and Starkman (2003), 
Saini et al. (2004)] that consider using the data to select the parametric 
model for w; we discuss these further in Section 6. 

This paper makes two contributions to this line of work. First, we intro- 
duce a new technique for testing hypotheses about w that does not require a 
specific parameterization of w. The technique is based on combining shape 
constraints on r, features of the functions in the null hypothesis, and any 
desired cosmological assumptions. As we show, this technique can be used to 
distinguish among currently competing models. Second, we develop a frame- 
work for nonparametric estimation of w with corresponding assessment of 
uncertainty. Given a sequence of parametric models for w of increasing di- 
mension, we use the forward operator T(-) to convert it to a sequence of 
models for r and use the data to select among them. We also show how to 
construct a representation of w that gives good performance in the forward 
mapping approach. 

Both of these methods take a fully nonparametric approach because, de- 
spite many ideas and some theoretical guidance, little is known about the 
function w. Although for current data low-dimensional models for w appear 
sufficient, the same may not be true with future data sets, which will have 
the precision to detect subtle structure. Accounting for model uncertainty 
in a rigorous and efficient way is one of the main values added by this work. 
Such uncertainty is not accounted for by the "figures of merit" proposed by 
cosmologists, for instance, in Albrecht et al. (2006), Sarkar et al. (2008), and 
Wang (2008), but it should be. For estimation within a single parametric 
model, we effectively use a maximum likelihood approach, so our inferences 
will have a precision comparable to the best methods in current use. But 
given that the structure of w is unknown and largely unconstrained a priori, 
some of that precision must be sacrificed to capture that structure at the 
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proper level of complexity. The methods we propose here also satisfy sev- 
eral other needs that cosmologists have. They can straightforwardly be used 
in combination SNe data with other types of data that probe dark energy. 
They allow great flexibility in the models and parameterizations used for w 
and in the cosmological assumptions that can be imposed. They provide an 
assessment of uncertainty on the inferences. And they are computationally 
efficient enough to handle the forthcoming large data sets. 

3. Data. Stellar magnitude is a logarithmic scale for the brightness of as- 
tronomical objects, defined so that dimmer objects have larger values. The 
apparent magnitude of an object describes how bright the object appears 
from Earth: m = 2.51og^Q(///o), where / is the flux of light produced by 
the object that is received by a detector on Earth (in some specified range of 
wavelengths) and /o is the corresponding flux for a reference object. The ab- 
solute magnitude of an object describes how bright the object appears from 
a fixed reference distance (10 parsecs), and thus is related to the apparent 
magnitude by M = m — 51ogiQ((iL/10), where di is a distance to the ob- 
ject measured in parsecs. [Specifically, is called the luminosity distance, 
which is one of several metrics used by astronomers; cf. Hogg (2000), and 
which differs from comoving distance by a factor of 1 + z.] The difference 
between them, /i = m — Af = 51ogio(dL/10), is called the distance modu- 
lus, a logarithmic measure of the distance to the object. For Type la SNe, 
apparent magnitudes are observed directly, and absolute magnitudes are de- 
termined from the observed luminosities over SNe lifetimes by fitting these 
observed "light curves" to a template and estimating their peak luminosities 
[Wood-Vasey et al. (2007)]. 

We analyze data for 192 SNe la from Davis et al. (2007) [see also Riess 
et al. (2006) and Wood-Vasey et al. (2007)]. The data includes, for each 
supernova: (i) redshift z, (ii) distance modulus /i, and (iii) standard error 



r = y + for the distance modulus, where is the intrinsic uncertainty 
in the distance modulus and is an estimate of error induced by the super- 
nova's peculiar velocity relative to a local standard of rest. (We ignore the 
uncertainties of the redshift estimates, which are generally less than 1%.) In 
the Supplementary Material [Genovese et al. (2009)], we also identify two 
supernovae whose data exhibit nontrivial influence in the analyses, but we 
include them because they could not be disqualified on concrete grounds. 

Let Ui and Zj denote the observed distance modulus and redshift, respec- 
tively, for the ith. supernova, i = 1, . . . ,n, where n = 192 and where, contrary 
to common astronomical practice, we have ordered that data by the redshift 
rather than date of supernova. We model Ui as Gaussian with mean n{zi), 
that is. 




(18) 
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where the e^s are assumed independent, mean zero, Gaussian noise terms 
with unit variance and the r^s are the given standard errors of the distance 
moduh measurements. 

We express the data in terms of comoving distance (assuming a flat uni- 
verse) by transforming as follows: 

(19) — ^ l0{U.-25)/5 ^ ^/^N . iQ{n/5)e, ^ 

C{l+Zi) 

where c is the speed of light. Thus, letting Yi denote the log^Q of the left-hand 
side of (19), we have 

(20) Yi = log-i^Q r (zi) + aiEi, i = l,...,n, 

where = rj/5. Figure 2 shows 10^ plotted against z with associated error 
bars; such a plot is called a "Hubble diagram." We thus call r "observable" 
because it can be directly estimated from the observed data. 

4. Nonlinear inverse problem formulation. We distinguish two uses of 
the word model in this paper. A cosmological model for dark energy is a set of 
assumptions about the underlying physics that gives rise to a particular form 
of the equation of state. A statistical model for it; is a family of probability 
distributions for the data indexed (at least) by a parameterization of w, 
possibly infinite dimensional. 

A particular cosmological model can be analyzed under a specific statis- 
tical model, but the scope of the inferences is limited by the viability and 
flexibility of the assumptions made. We consider statistical models whose 
stochastic component is specified by equation (20); each such model is then 
determined by the parameters Hq and Qm and a representation of w. 

We now re-express the relationship between the comoving distance and 
the equation of state as an explicit analytic expression mapping w (and 
the cosmological parameters Hq and 0^) to r. Equation (11) describes the 
relationship between comoving distance and redshift. By expanding H(z) as 
in Huterer and Turner (2001), one can derive the following equation: (we 
give an alternative derivation in the on-line Supplementary Material) 

r{z) = H^^ r ds[nrr,{l + sf + (1 - nm){l + sf 

(21) 

X e"^ /o -">(«)/(!+«) d«j -1/2 _ 

We can thus write r as the image of an operator acting on w; specifically, 
we can write r = 10^('"'-'^0'^'"), where the nonlinear operator T{w; HQ,0,rn) 
is defined by logj^g of the right-hand side of equation (21). Equation (20) 
thus describes a nonlinear inverse problem 

(22) Yi = T{w;Ho,nm){zi) + aiEi, i = l,...,n. 
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5 




Redshift z 

Fig. 2. The data from equation (20) transformed to plot comoving distance versus red- 
shift. Standard error bars are attached to each point. 



The operator T is complicated and depends on two unknown parameters, 
but it does have several useful properties, as we will see below. An important 
use of equation (21) is to translate a model for the unobservable w into the 
observable r. For instance, given any parameterization of w, equation (22) 
determines a likelihood function. We show how to use this for parametric or 
nonparametric inferences in the following sections. 

It is also sometimes helpful to consider the universe's acceleration di- 
rectly. Cosmologists traditionally express this in terms of the dimensionless 
deceleration parameter^^ q{^)^ which is defined in terms of the scale factor 
by 



(23) 



qit) 



d{t)a{t) 



-1 



HHty 



^This was named when it was thought that the universe was decelerating. 
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and in terms of redshift and comoving distance by 

(24) g(z) = -l-(l + z)^. 

Using a similar method to the above, we can solve for q in terms of w and 
for r in terms of q to obtain 

(25) q{z) = \ + ?,w{z) , 

(26) r{z) = Fq-i r ds e- ^ duii^ii^)) I (x+u) _ 



These equations are derived in the Supplementary On-line Material. 

Equations (21), (25) and (26) have several properties that are valuable for 
statistical inference. First, note that for any tt;, r'(z) > 0, so r is a monotone 
increasing function of z with r(0) = 0. In fact, 

(27) ii+l <r'(z)<ii±^, 

where the upper bound can be made sharper if w is assumed bounded from 
below. Second, r is monotone decreasing in w for each fixed value of 
and ^rn- Specifically, if w\ and are two candidate equations of state 
with corresponding comoving distance functions r\ and r2, and if 1^2 (-z) ^ 
tL'i(z) for all z > 0, then rii^z^ < ri{z) for all z>0. Similarly, q is monotone 
increasing in w; W2{z) > wi{z) for all z > implies that q2iz) > qi{z) for 
all z >0. And r is monotone decreasing in q; q2{z) > qi{z) for all z > 
implies that r2{z) < ri{z) for all z>0. Third, as shown in the Appendix, for 
r to be concave it is sufficient that w{z) > —1/(1 — ^m) for all z>0. Under 
mild smoothness assumptions on w, the concavity of r holds more broadly. 
Fourth, in both equations involving w, taking i(; = is equivalent to taking 

For any specific parameterization of w and any choice of and Hq, it 
is straightforward to evaluate r numerically. For instance, under a constant 
w model w = wq, equation (21) reduces to 

(28) r{z) = H^' r dsl^rail + sf + (1 - nm){l + s)3(l+-o)]-l/2. 

JO 

More generally, expanding w{z) = —J2j Pji^ji^) ^ (not-necessarily or- 
thonormal) basis -00 , V'l ; • ■ • yields 

(29) r{z) = H^' rd5[o^(i + 5)3 + (i_j7^)(l + 5)3g-3E,ft^.W]-i/2^ 

Jo 

where ^j(s) = Jq ipj{u)/{l + u)du. Taking the expansion to be finite gives 
three important special cases: 
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1. Polynomial in z: il)j{z) = z^ , j = 0, . . . , d, giving 

r{z)=H^^ rds[nm{l + sf 

(30) 

+ (1 _ a„)(l + >,)3(l-ao)e-3E,.i(-l)^°.-Vipl/2^ 

where ak = Y.j=ki-^y Pj for A; = 0, . . . , d. 

2. Polynomial in the scale factor a: ipj{z) = (1 + z)~^ , j = 0, . . . ,d, giving 

r{z)=H^^ rds[n^{l + sf 

(31) 

+ (1 - + ,)3(l-/5o),3E,../3,((l+.)-.'-l)/,pl/2_ 

3. Piecewise constant: ipj{z) = ^{sj,sj_^-^]iz) for j = 0, . . . ,K — 1, where = 
So < si < • • • < s;^ are breakpoints for K fixed bins and where '^-(sj^Sj^-^] (z) 
is 1 if Sj < z < Sj+i and otherwise. In this case, equation (21) becomes 



■3B(s)]-l/2 



(32) r{z) = H^^ r ds[^lm{l + sf + {l-VL^){l + sfe 

JQ 

where 

J(s) 

(33) i?(,) = ^/?^.log(-i±i^) +/J iog(,A±^ 

and where J(s) = max{0 <j<K: Sj < s}. Despite the discontinuities in 
w, this expression is a smooth function of the /3 parameters. 

Extension to other bases — such as B-splines, orthogonal polynomials, and 
wavelets — is straightforward. 

Combined with equation (20), each of these expressions produces a like- 
lihood for w, Hq, and f^m- Although nonlinear, these likelihoods are well- 
behaved for optimization purposes, and weighted, nonlinear least-squares 
is computationally efficient in practice. Good estimates of the coefficients 
can be obtained for a wide variety of models, which in turn supports both 
parametric and nonparametric inferences about w. 

5. Methods I: hypothesis testing. As discussed above, one method for 
distinguishing among models of dark energy is to first estimate the equation 
of state and use this estimate to test hypotheses about cosmological models 
[e.g., Huterer and Cooray (2005)]. This approach has several disadvantages, 
including that the power of the test depends on having a good estimator and 
that it requires accurate standard errors for the entire function. Moreover, in 
practice, such tests usually rely for their validity on an assumed parametric 
model for w. It would be desirable to be able to test cosmological models 
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without a preliminary estimator or assumed parameterization, and in this 
section, we construct a method to do that for certain classes of hypotheses. 

The basic idea is that we use the forward operator given in equation (21) 
to map a set of possible w's to the r domain, and use the data to test the hy- 
pothesis there by inverting a nonparametric confidence set for r. Two issues 
arise in such a scheme. First, for general sets in the w-domain, it can be diffi- 
cult to compute their image in the r-domain, but we use the properties of the 
operator discussed in Section 4 to easily compute the mappings for certain 
classes of hypotheses. Second, performing a sharp nonparametric test (or 
constructing a small nonparametric confidence set) can be difficult without 
structural assumptions, but we take advantage of strong shape constraints 
satisfied by the comoving distance function. 

Here, we consider null hypotheses of the following forms: 

A. simple equalities for w. w = wq, 

B. inequalities ioi w. wq < w < wi, 

C. inequalities for w': Wq<w'< w'l, 

D. inclusion: w £V for a linear space V of fixed dimension, 

and various intersections of these, where wq, wi, w'q, and w'l denote various 
fixed functions, not necessarily constant. [We use the inequality w <wo to 
mean that w{z) < wq{z) for all z, and similarly for other inequalities between 
functions.] 

Testing such hypotheses gives direct tests of various cosmological models. 
The null hypothesis that the cosmological constant model holds, for example, 
translates to a simple null hypothesis with wq = —1. Quintessence solutions 
lead to a variety of constraints on w and w' that can be tested by combining 
hypotheses that are inequalities for w and for w' . For instance, as we show 
in the Appendix, thawing solutions satisfy 

(34) l+y;<^<3(l+y;), 

a ma 

and freezing solutions satisfy 

(35) 3w(l + w)<-^<0.2w(l + w), 

a In a 

when — 1 < u;<-0.8 [Caldwell and hinder (2005)], where a is the scale factor. 
These bounds can be re-expressed for w in the same range as 

m i±^-i, 

{l + z}-^ 1 + z 

and 

'u;(0) 



(l + z)3+«;(0)((l+z)3-l) 

(37) 

< w[z) < 



(l + z)0-2+u;(0)((l + z)0-2- 1)' 
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where w{0) is a free parameter. 

The strategy underlying our testing procedure is to use equation (21) to 
translate hypotheses about w into hypotheses about r, making it possible 
to test any hypothesis that translates into a manageable form. Our test is 
derived by inverting a 1 — a confidence set for {r{zi), . . . ,r(z„)), as follows: 

0. Select a small < a < 1. 

1. Construct a 1 — a confidence set C for the unknown vector {r{zi), . . . , t(z„)) 

2. Construct the set Rq of vectors (ro(zi), . . . , ro(zn)) where tq is a co- 
moving distance function produced by an equation of state consistent 
with the null hypothesis. 

3. Reject the null hypothesis if C H i?o = 0. 

In practice, the sets in steps 1 and 2 need not be constructed explicitly. For 
example, C typically takes the form of bands — an interval at each Zi — or a 
ball centered around a particular vector. And Rq can usually be represented 
implicitly by an efficient search over the null hypothesis in w (mapped for- 
ward by T) targeting those r's that lie outside C. In practice, this procedure 
can be made computationally efficient for a broad range of hypotheses. 

One way to define the confidence set C is the set of vectors for which a 
standard chi-squared goodness-of-fit test does not reject the null hypothesis. 
In light of equation (20), the chi-squared goodness-of-fit ball gives a confi- 
dence set for {logio{r{zi)), . . . ,logiQ{r (zn))), which is easily transformed into 
a confidence set for (r(zi)), . . . ,r(z„)). As the number of data grows, how- 
ever, the chi-squared confidence sets become unduly conservative, reducing 
the power of the test. So, we also use alternative confidence set procedures 
that produce smaller confidence sets, giving the test higher power [Baraud 
(2004), Davies et al. (2007), Ingster and Suslina (2006)]. 

Suppose that for i = 1, . . . ,n, Yi = fi + aiSi, where the <TjS are known 
numbers and the e^s are independent Gaussian variables. This corresponds 
to equation (20) with fi = logi(jr{zi). If the vector / = (/i, . . . , /„) denotes 
the true but unknown values of the function at the observed points, then a 
1 — a confidence set C for / is a random set, constructed from the data, that 
satisfies 

(38) P{C3f}>l- a. 

We want C to be as small as possible. (Although there are several reasonable 
definitions of size, we will use the simplest: the radius of a confidence ball 
in the corresponding norm, the width of confidence bands, and the volume 
of a more general set. Our methods work as well if the confidence set is 
constructed to optimize some other criterion.) 

There are several ways to construct such confidence sets. One way is to 
invert a chi-squared goodness-of-fit test for the null hypothesis f = f^, giving 
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C = {f:T\fo)<xlJn}, where 

Xn.a is upper-tail a quantile of the 
corresponding chi-squared distribution and where r^(/o) = (l/'n,) X^iLil^ ~ 



The chi-squared confidence set is simple to use, but it has several ma- 
jor drawbacks. The confidence set is relatively large; the radius of the set 
Xn,a/Vn is 0(1) no matter how large n is. The set is constructed from a 
rough estimator of /, namely, the data. The size of the set is independent 
of the data and thus cannot adjust to evidence of smoothness. And some 
prior information, such as shape restrictions, is difficult to incorporate in 
practice. 

There are practical confidence set procedures that address all these draw- 
backs. We consider two: shape-restricted confidence bands from Davies, Ko- 
vac and Meise (2007) and adaptive chi-squared confidence sets from [Baraud 
(2004)]. Both provide finite-sample, honest [Li (1989)] confidence sets that 
adapt their size based on the data. Both provide substantially smaller con- 
fidence sets than the chi-squared ball. Both are computationally practical, 
though somewhat more work than the naive chi-square confidence set. And 
both allow us to incorporate prior information about the comoving distance 
to produce a smaller confidence set. Because the Davies et al. procedure 
performed better in simulations with current sample sizes and because it 
requires fewer tuning parameters, we focus on that procedure in this paper. 
The shape constraints are also well adapted to our prior information about r. 
In principle, the Baraud procedure combined with shape constraints should 
outperform the Davies et al. procedure, but we will explore this comparison 
in a future work. More detail about the Baraud procedure is given in the 
Supplementary On-line Material. 

The procedure of Davies et al. generates confidence bands under the as- 
sumption that / is monotone and concave. For each 1 <i <k <n, define the 
integer interval Jj^ to be the set of indices 1 <j <n such that Zi < Zj < z^. 
Thus, if there are no ties among the redshifts, lik = {i, . . . ,k}. Consider 
the following statistics, modified from Davies, Kovac and Meise (2007), to 
account for the different standard errors of the measurements: 



where ^(•) gives the cardinality of a set. These serve as test statistics for 
testing whether the residual mean is zero along any index interval. When 
/ is the true vector, the Tjfc(/)s each are mean zero Gaussian variables. 
Davies, Kovac and Meise (2007) point to a procedure for computing these 
statistics in O(nlogn) time and offer an approximating subset of size 0{n) 
for large n. We use the latter for convenience, but use of the full set did not 
significantly change the results. 




(39) 
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The procedure begins with a confidence set for the Tj^'s. In Davies, Ko- 
vac and Meise (2007), this is a confidence cube with edge length equal to 
twice the 1 — a quantile of maxj ,t |Tjfc(/)| . The key for constructing the final 
confidence set is that the initial confidence set has linear boundaries, mak- 
ing consistency with the initial confidence set a computationally tractable 
constraint. Constraints for concavity and monotonicity are also linear. This 
gives a confidence set for / consisting of those vectors g such that the Tjk{g)s 
lie in the initial confidence set and that g satisfies the imposed shape con- 
straints. We could use this confidence set directly, but it is computationally 
much simpler with the above types of hypotheses to use confidence bands. 
We compute the confidence bands are constructed by maximizing and min- 
imizing fj subject to the vector / lying in the initial confidence set and 
satisfying the shape restrictions. These require optimization of a linear ob- 
ject function with linear constraints and thus can be solved with two linear 
programs (i.e., optimization problems with a linear objective and linear con- 
straints) for each j. 

We modify the Davies, Kovac and Meise (2007) procedure in several ways. 
First, because monotonicity and concavity are used in the procedure in log 
space, the confidence bands need not be concave in r space. We adjust for 
this by optimizing the bands in r space, finding the smallest bands consistent 
with the shape restrictions there, including the additional constraint that 
r(0) = 0. This involves two additional linear programs for each Zi. Second, 
we can use a smaller initial confidence set for the Tj^s with some additional 
computation. Specifically, the distribution of Tjfc(/)'s is a degenerate Normal 
whose covariance depends only on the collection of /jfc intervals. With an 
eigen-decomposition of this covariance, we replace the hyper-cube of Davies, 
Kovac and Meise (2007) by a substantially smaller degenerate ellipsoid and 
get proper coverage. The 0{n) sized symmetric eigen-decomposition is ex- 
pensive for large n but can be parallelized if necessary. This also requires 
optimizing a linear function subject to linear and ellipsoidal constraints. 
Such a problem can be reduced to a convex optimization problem called 
second-order cone programming [Boyd and Vandenberghe (2004)], in which 
we minimize a linear function subject to one or more "second-order 
cone constraints" of the form \\Ax + 6|| < c^x + d, where || • || is the Eu- 
clidean norm, ^ is a matrix, b and c are vectors, and d is a scalar. Many 
common convex optimization problems can be reduced to this form. We 
used the MOSEK and CPLEX software^^ to compute the solutions to these 
problems. Finally, to test whether a linear space of vectors intersects the el- 
lipsoidal constraint set and satisfies the shape constraints, we can minimize 
the (quadratic) distance to the center of the ellipsoid subject to the linear 



See http://www.mosek.com and http://www.ilog.com/products/cplex/. 
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shape constraints and inclusion in the given space. If the resulting optimum 
is sufficiently close to the center, the test is affirmative. This is yet another 
form of convex optimization (called quadratic programming) that can be 
implemented with the same software. 

In all of the hypotheses we test, there are one to three free parameters 
that must be varied, usually including Hq and Qm- A simple grid search 
is practical and straightforward in these cases. For any value of the free 
parameters, we can tell whether the corresponding r lies in the confidence 
bands by direct comparison. We can also incorporate information from other 
studies by using confidence sets on {Hq, Qrn) to define this subsidiary search. 
Assuming independent data sets and taking a' = 1 — \/l — a, if we use a 
1 — a' confidence set for the cosmological parameters and in our procedure, 
the resulting test has level a as required. 

Our procedure works also under more restrictive assumptions about the 
form of w, with correspondingly sharper results as the assumptions grow 
stronger. For this, the confidence set C is constructed using the assumed 
parameterization. The resulting test will have higher power than the non- 
parametric test when the assumed parameterization holds. Note, however, 
that the validity of any inferences under a specific parameterization depends 
strongly on the parameterization being accurate. For instance, we assume 
here that w > —1 [Carroll, Hoffman and Trodden (2003)], which implies that 
r is monotone concave. But the same basic procedure works, with somewhat 
lower power, when that assumption is dropped. 

Step 2 of the procedure depends specifically on the hypothesis being 
tested. We now derive the sets i?o for null hypotheses of the forms listed 
above. Let M. denote the set of vectors {r{zi), . . . ,r{zn)) for functions r 
that meet the a priori conditions that the comoving distance must satisfy: 

A. Under a simple null hypothesis w = wq, equation (21) generates a two- 
parameter family of functions tq as Hq and vary; Rq is the set of 
vectors (ro(zi), . . . ,ro{zn)) for tq in this family. 

B. Under the null hypothesis, w > wq, equation (21) shows that, for fixed 
Hq and ^Im, r < vq, where vq is produced in (21) hy w = wq for the given 
value of Hq and 0^. Again, varying Hq and Qm produces a two-parameter 
family of functions vq. Rq is the set of vectors (ri, . . . , r„) G Ai such that 
fi ^ ^oizi), . . . , rn> i^oi^n) for some vq in the family. The restriction to 
A4 sharpens the results. It is not strictly necessary, but because equation 
(21) produces functions in M, it is an improvement that is virtually 
cost free. The other direction of inequality is handled similarly, using the 
monotonicity of r in 

C. Null hypotheses of the form w' > w'q can be handled by re-expressing 
the exponent in equation (21). Integrating by parts and writing 'w{s) = 
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w{0) + Jq w'{u) du yields that 

(40) / -^^du = w(0)log(l + s)+ / w'(u)(log(l + s) -log(l + u))du. 
Jo l + u Jo 

The second term in the right-hand side integrand is nonnegative, so w' > 
w'q imphes, for fixed Hq, and w{0), that r < tq, where ro is the right- 
hand side of equation (21) corresponding to {'w{0),Wq, HQ,Qrn)- Varying 
Hq, Qjji, and w{0) produces a three-parameter family of functions, and 
as before, Ro is the set of vectors in M whose components are at least 
as big everywhere as some function in this family. Other inequalities in 
w' are handled similarly. 
D. The null hypothesis that w lies in some linear space of functions V is 
useful primarily to test the goodness of fit of statistical models for w. 
We select an arbitrary basis for V and form a dim(l/) -|- 2 dimensional 
family of functions corresponding to each (ifoj^m) and each vector of 
coefficients in the basis expansion. Rq is the set of vectors produced by 
these functions evaluated at zi,...,Zn- See equation (29). This case is 
handled in practice by numerical optimization and thus works best for 
low to moderate dimensional spaces. It is not necessary to restrict to a 
linear space, but that is the best behaved case numerically. 

Note that this same approach can be used to test hypotheses about the 
deceleration function q{z). For instance, we may wish to test the null hypoth- 
esis that the universe is nonaccelerating. This can be expressed in several 
ways. First, we can test whether the universe is matter dominated with- 
out dark energy, which corresponds to w = 0, or equivalently, = 1, or 
q = 1/2 from equation (25). We call this the matter dominated hypothesis. 
In contrast, a nonaccelerating universe corresponds to q > 0; we call this 
the strongly nonaccelerating hypothesis. Using the monotonicity properties 
of equation (26), this matter-dominated hypothesis maps to the null hy- 
pothesis r = Hq^{2 - 2(1 + z)"^/^) and the pure nonaccelerating hypothesis 
corresponds to the null hypothesis r > H^^ log(l -|- z). Both generate a one- 
parameter family and corresponding Rq. Note two issues with the latter, 
one-sided hypothesis. By taking Hq large enough, we can make the lower 
bound as small as possible and the null trivially true, so we need to use a 
prior confidence set for Hq as described above. Second, the condition q>0 
allows expansion histories that are strongly at odds with current theory. So, 
strictly speaking, testing q>0 will offer poor power against alternatives we 
care about. We fix this problem in two ways: combining q >0 with (i) a 
bound on w such as w > —1/3 to preclude acceleration in a dark-energy 
dominated universe, and (ii) assume that q exhibits a change point between 
the matter-dominated case q = 1/2 and the dark-energy dominated case 
q >0. Both tighten the bounds substantially, at the expense of an added 
restriction or free parameter. 
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6. Methods II: estimating the equation of state. Answering the main 
questions about dark energy that are currently being addressed with su- 
pernova data involves testing among competing models. But ultimately we 
want to estimate the equation of state. As described earlier, this problem 
has received substantial attention in the literature. In this section we de- 
scribe a framework for constructing nonparametric estimators of w with an 
associated assessment of uncertainty. 

A nonparametric procedure for estimating the equation of state treats 
w as an infinite parameter belonging to a specified space of functions. Of 
course, with only a finite amount of data, any estimator has limited resolu- 
tion, so we must use the data to determine the complexity of the estimator. 
This gives rise to the bias-variance trade-off: too high a complexity pro- 
vides a better apparent fit but gives estimates with high variance; too low a 
complexity gives estimates with low variance but bias from model misspec- 
ification. 

There have been several related works in the dark energy literature. Saini, 
Weller and Bridle (2004) formulate the problem as a Bayesian inference prob- 
lem where the parameter space includes disjoint spaces of low-dimensional 
polynomials. Huterer and Starkman (2003) construct an empirical basis via 
principal components. They correctly note that choosing the number of ba- 
sis elements in a nonparametric analysis requires that one carefully balance 
the bias and the variance. We provide a concrete mechanism for choosing 
the number of elements in practice. 

Our approach begins with a collection of models M.i,M.2, ■ ■ ■ , where each 
Alfc is a linear space of functions with dimension k. We then select a model 
that balances bias and variance by minimizing an empirical measure of risk 
or comparable criterion. We use BIG [Bayesian Information Criterion, also 
called the Schwarz criterion; Schwarz (1978)], as it is both simple and ef- 
fective, but other criteria give similar results. All the spaces Mk — including 
those for small k — contain smooth and constant functions, and as k in- 
creases, the spaces add detail to capture more complex fluctuations. We 
define J^k to be the fe-dimensional space of cubic B-splines [de Boor (2001)] 
over the range of the data, with equally spaced knots. The choice of knots 
could be optimized for even better performance by adapting it to the dis- 
tribution of redshifts in the data. Note that with this choice, the bases for 
the Aik spaces are not nested. 

A common alternative choice in nonparametric function estimation is to 
choose an orthonormal basis to represent w and define Mk to be the span of 
the first k terms. In an inverse problem such as this, however, this strategy 
requires care. Choosing a basis to obtain an efficient estimator requires bal- 
ancing the information passed by the forward operator and the conciseness 
of the representation for w. For example, the Fisher basis [e.g., Huterer and 
Cooray (2005)] based on the forward operator in equation (22) has most of 
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Fig. 3. The first four Fisher basis functions, all of which decay quickly to zero in redshift. 



its variation at low to moderate redshift where w is relatively more informa- 
tive in the data. But as shown in Figure 3, the basis elements decay quickly 
to zero. To fit a simple smooth function, even a constant, requires a large 
number of basis functions, leading to high variance estimators. This is an 
especially important issue when considering variations around the cosmo- 
logical constant model. 

It was for this reason that we defined the Mk as above. Each con- 
tains constant and smooth functions over the entire redshift domain. As 
a consequence, a good estimate of a smooth function can be obtained in 
any of the models, and when w is in fact smooth, the procedure will select 
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a lower-dimensional model, with consequent gains in mean squared error. 
Moreover, unlike polynomials, the B-splines are localized and numerically 
well conditioned. 

Of course, for most purposes, an estimate of w is not sufficient, so we 
also need an assessment of uncertainty in the estimates. Consider first in- 
ference for a single Mk- It is possible to use a Normal approximation at the 
maximimum likelihood estimate from the nonlinear regression to derive an 
approximate covariance matrix for the estimated parameters. In practice, 
this seems to perform well, but because it is difficult to bound the accuracy 
of the approximation given the nonlinearity, we use a resampling approach. 
The error bars can be computed using a parametric or nonparametric boot- 
strap. We prefer to use a nonparametric bootstrap since the error bars will 
then be less affected by any bias in the specification of the model, especially 
given that we are selecting among models initially. We generate bootstrap 
confidence intervals for the parameters in the model by resampling residuals 
from the model fit, renormalized to have appropriate variance [Efron (1979), 
Efron and Tibshirani (1994)]. The basic procedure is as follows: 

1. Compute the maximum likelihood estimator {(3,9), where /? is the vector 
of parameters for w and = {Hq,^}^)- 

2. Compute residuals Cj = 1^ — f{zi; f3, 9). 

3. Using the standard errors of Yi and the linear approximation at the maxi- 
mum likelihood estimator, standardize the residuals to unit variance. Call 
these standardized residuals Si. 

4. For h = l,. . . ,B, for some large B, draw pseudo-noise from the empirical 
distribution of the Sj. Call these for i = 1, . . . , n. 

5. Generate pseudo-data 



7. Compute standard errors and confidence intervals for these parameters 



from the (/?*(*), ^*(''))'s as in Efron (1979) and Efron and Tibshirani 



We use the bootstrap confidence intervals to compute confidence bands for w 
and q by computing the largest and smallest values of the functions at each 
redshift that are consistent with the confidence intervals on the parameters. 

In the nonparametric case, it is common practice to use the confidence 
bands corresponding to the selected model. These are straightforward and 
accurate when the selected model holds but are necessarily optimistic be- 
cause the bands do not account for the variation in the model selection 
process or for the potential bias induced by choosing too simple a model. 




(1994). 
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A simple improvement we use is to incorporate the model selection into 
the resampling process, using the largest and smallest estimated w or q 
from the bootstrap samples to construct the confidence bands. Because this 
effectively includes model bias in the bands, this approach is likely to be 
somewhat conservative. We will explore other approaches to this problem in 
a future paper. 

7. Results. 

7.1. Testing cosmological models. We test seven cosmological models us- 
ing the procedure described earlier, independently of any parameterization 
for w. Three models (cosmological constant, frustrated cosmic strings, and 
domain walls) can be tested with a simple null hypothesis of the form w = wq. 
The two quintessence models (thawing and freezing solutions, resp.) were 
tested with inequality null hypotheses given by equations (36) and (37) in- 
tersected with the condition that —l<w< —0.8. [We also tested more ex- 
pansive versions of these hypotheses using (i) equations (36) and (37) alone 
and (ii) the hypothesis if' > and if' < intersected with the condition that 
— l<w< —0.8. But being strict supersets of the original null hypotheses, 
these are less likely to reject.] We tested both the matter dominated and 
strongly nonaccelerating universe hypotheses. For the latter, we used a con- 
fidence interval obtained from the current best estimates^^ and adjusted the 
confidence level as described in Section 5. Finally, we tested the inclusion 
hypothesis that w is a constant, possibly different from —1. Table 1 shows 
the results of these tests at various significance levels. The no dark energy 
model is clearly inconsistent with the data (p-value p~ 0), but none of the 
other models are rejected at the 13% level. Note, in particular, that the 
cosmological constant is consistent with the data. 

A false null hypotheses might fail to be rejected because the power of the 
test is too low. Because our procedure has essentially as much power as possi- 
ble given the available information about the only ways to improve power 
are either to make stronger assumptions about the form of w or to get more 
data. We argue that the latter is necessary. The results do not change when 
performing the same tests assuming a linear form w{z) = —{Po + (^iz), which 
is the simplest nontrivial parameterization and a correspondingly smaller 



^ For current estimates of cosmological parameters, see 
http://lambda.gsfc.nasa.gov/product/map/current/. We adopt the value of H^Q.^ 
derived from "all" data for the LCDM model. This gives Ho 72 ± 8 km/s/Mpc. Note, 
however, that the absolute value of the Hubble constant is only determined by the 
supernova data up to an arbitrary shift because of calibration of the absolute magnitudes. 
We thus recentered this confidence interval around 65, which is consistent with the 
calibration of our data. 
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Table 1 

Results of nonparametric hypothesis tests for various cosmological models of w. The 
significance levels correspond to 1, 1.5, 2, and 2.5 standard deviations respectively from a 

Gaussian mean 



Model 




Rejected at level 




32% 


13% 


5% 


1% 


Cosmological constant 


yes 


no 


no 


no 


Frustrated cosmic strings 


yes 


yes 


yes 


no 


Domain walls 


yes 


no 


no 


no 


Matter dominated 


yes 


yes 


yes 


yes 


Nonaccelerating 


yes 


yes 


yes 


no 


Quintessence thawing 


no 


no 


no 


no 


Quintessence freezing 


no 


no 


no 


no 


Constant w 


no 


no 


no 


no 



confidence set. The pattern of rejections is basically the same, and in par- 
ticular, there is insufficient evidence to move away from a cosmological con- 
stant. Of course, there is no reason to believe the linear form for w, and if it 
is false, inferences under that assumption can be misleading. But this shows 
that strengthening the assumptions is not enough to overcome the lack of 
information in the data. 

7.2. Fitted models for the equation of state. Our fitted w is a constant, 
w = -1.013 ± 0.124 with = 0.268 ± 0.028 and Hq = 65.6 ± 0.90, where 
the standard errors are based on 1000 bootstrap iterations. (Note that the 
supernova data determine the value of the Hubble constant plus an arbitrary 
shift induced by calibration of the supernova absolute magnitudes, so it is 
the relative uncertainty rather than the absolute value that matters here.) 
The bootstrap BCa 95% confidence intervals do not differ much from the 
Normal intervals based on the bootstrap standard errors: [—1.262,-0.796] 
for w, [0.220,0.324] for 0^, and [63.9,67.4] for Hq. It is straightforward to 
compute joint confidence sets for these parameters, but we do not report 
them here. As a check on these results, we note that within the parametric 
models we consider (polynomials in z or a, piecewise constants, B-splines) 
likelihood ratio tests between the constant model and the higher-order mod- 
els in the family fails to reject with p-value p > 0.85. In all three cases, BIG is 
monotone increasing with the constant model the clear choice. A likelihood 
ratio test of the cosmological constant versus the constant w model fails to 
reject with j*- value p = 0.137. 

7.3. The need for more data. A key question is whether current super- 
nova data are sufficient to resolve the differences among interesting models 
for w and q. We argue here that the answer is no. 
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Table 2 

Power for distinguishing a constant model w = wi from a cosmological constant w = —1 
using a likelihood ratio test with significance level a. Each is based on 3000 simulations 

of data from equation (20) 



Alternative (lOi) 




Significance 


level (a) 




0.68 


0.87 


0.95 


0.99 


-0.85 


0.352 


0.169 


0.080 


0.015 


-0.90 


0.328 


0.157 


0.056 


0.012 


-0.95 


0.303 


0.121 


0.047 


0.012 


-0.99 


0.322 


0.116 


0.044 


0.007 



First, even under strong assumptions and with essentially optimal proce- 
dures, there is not enough evidence to distinguish among interesting models. 
The cosmological constant model is suggestively on the boundary at the 13% 
level, but no conclusive differences are supported by the data. 

Second, we can use as a minimal criterion for resolvability the power of the 
likelihood ratio test for distinguishing the cosmological constant model from 
a constant w model. With existing standard errors, it is straightforward to 
compute this power by simulation through equations (21) and (20). Table 2 
shows the power of this test for various significance levels and alternatives, 
all of which are low. The power for distinguishing a constant w model from a 
piecewise constant model with one breakpoint are lower for a similar variety 
of alternatives. 

Third, an even more striking demonstration of model degeneracy is given 
by Figure 4. This shows two very different equations of state that give virtu- 
ally indistinguishable fits to the data, with a chi-squared deviation of 0.04. 
This example is driven primarily by uncertainty in $7^ , which can be reduced 
using other (non SNe) data. 

8. Discussion. Matching the many studies in the astrophysics and cos- 
mology literature, we find that current supernova data do not yet provide 
tight enough inferences to make nontrivial claims about w. And in particu- 
lar, the data are consistent with a cosmological constant model. This puts 
the focus squarely on future data, which will likely provide notably stronger 
constraints on the dark energy equation of state. This raises the question of 
what this newfound precision will mean. If the cosmological constant model 
holds or if w is very close to but not equal to —1, then we will be in much the 
same state as we are currently. But otherwise, the insight we gain into the 
nature of dark energy will depend on our ability to distinguish competing 
models and to infer subtle structure in w. For this purpose, a nonparametric 
approach will be particularly effective. This paper describes a new technique 
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CO 




0.0 



1.0 



1.5 



Fig. 4. Two equations of state that produce statistically indistinguishable r's. 

for testing among competing dark energy models with minimal assumptions 
and describes a framework for nonparametric estimation of w. 

Several challenges and open questions remain. Computation of sharp, hon- 
est confidence sets in nonlinear inverse problems is a mathematically diffi- 
cult problem but important to getting as much information as possible from 
the data. Expanding the scope of the hypothesis testing methods to new 
classes of hypotheses that eliminate "unphysical" possibilities will improve 
the power of the technique. And systematic errors will become a major issue 
as the statistical uncertainties become smaller. 

The next generation of supernova data sets may answer many questions 
with the methods presented herein, but they will also raise a host of inter- 
esting new statistical problems. 



To derive equations (36) and (37), we begin by transforming equations 
(34) and (35) from the scale factor a to redshift z. Replacing dw/d\na by 
— (1 + z)w'{z) and reversing the inequality because of the negative sign yields 
the corresponding equations 



APPENDIX: EQUATION OF STATE BOUNDS FOR THE 
QUINTESSENCE MODELS 



(42) 




for thawing solutions, and 



(43) 



^.2w{z) 



l + w{z) 
l + z 



< w'{z) < —3w{z) 



1 + ^(2;) 
l + z 
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for freezing solutions. 

Begin with the assumption that w > —1, which we will weaken below. 
For the thawing equalities, divide through by 1 + if to get w' /{I + w) = 
(log(l +w)y and, thus, 

(44) - 3 J— < (iog{l + w{z))y < ^ 



1 + Z ~ ' ' V /// - -^^j^^ 

Integrating through from to z yields 

(45) - 31og(l + z)< log(l + w{z)) - log(l + w{0)) < -log(l + z), 
and taking exponents, 

(46) (l + z)"^<i±^<(l + z)-\ 

1 + w{0) 

which leads directly to equation (36). 

Similarly for the freezing solutions, w'/w{l +w) = (log(^i'/(l + w))y. Di- 
viding through and integrating as before gives 

(47) - 0.21„g(l + .) < log(^) - lo<r^) < -31og(l + 

Taking exponents and simplifying gives equation (37). 

For freezing solutions, w' > 0, and if w{0) = —1, then w = —1. For thawing 
solutions, either w{0) = —1 and w = —1, or w{z) = —1 for some z > 0. The 
latter case leads to a contradiction given the bounds on w' and continuity 
of w. Hence, the bounds hold for —l<w< 0.8. 

Acknowledgment. The authors would like to thank the referee for an 
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SUPPLEMENTARY MATERIAL 

On-line supplementary material for Genovese et al. 2009 

(DOI: 10.1214/08-AOAS229SUPP; .pdf). We provide further technical de- 
tails on data issues, derivations, and methods. 
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